#ifndef DYN_LYAPUNOVTIMESERIES_H
#define DYN_LYAPUNOVTIMESERIES_H

#include <armadillo>

using namespace arma;

/*! \author Joey Dumont <joey.dumont@gmail.com>
 *
 * \date 2012-08-24
 *
 * \brief Computes the Jacobian matrix and Lyapunov exponents of maps and
 * flows using a time series.
 *
 * Using only a time series has its input, this class computes the
 * the jacobian matrix at either a specified point of the orbit,
 * or computes the global Lyapunov exponents of the orbit.
 *
 * The approximation of the Jacobian is done via
 * a least-squares fit of an evolution operator.
 * After that, we use the formula given in Sano & Sawada (1985)
 * to compute the Lyapunov spectrum. We use the singular-value
 * decomposition to compute the orthogonal basis vectors.
 */

class LyapunovTimeSeries
{
public:
    /*! The constructor takes the time series and
     * the order of the Jacobian (we compute \f$J^p\f$).
     */
    LyapunovTimeSeries(mat _timeSeries, int _p, double _timeStep, double _epsilon, int _maxNumberOfPoints);

    /*! @name Accessor Functions
     * We access the variables of that class.
     */
    //@{
    mat getTimeSeries(){return timeSeries;}
    mat getJacobian(){return jacobian;}
    colvec getLyapunov(){return lyapunov;}
    int getP(){return p;}
    int getMaxNumberOfPoints(){return maxNumberOfPoints;}
    double getTimeStep(){return timeStep;}
    double getEpsilon(){return epsilon;}

    void setTimeSeries(mat _timeSeries){timeSeries=_timeSeries;}
    void setJacobian(mat _jacobian){jacobian=_jacobian;}
    void setLyapunov(colvec _lyapunov){lyapunov=_lyapunov;}
    void setP(int _p){p=_p;}
    void setMaxNumberOfPoints(int _maxNumberOfPoints){maxNumberOfPoints=_maxNumberOfPoints;}
    void setTimeStep(double _timeStep){timeStep=_timeStep;}
    void setEpsilon(double _epsilon){epsilon=_epsilon;}
    //@}

    /*! @name Main Methods. */
    //@{
    mat computeJacobian(int orbitPoint, int order);
    colvec computeLyapunovSpectrum(int order);
    //@}

protected:
    /*! @name Variables */
    //@{
    mat timeSeries;
    mat jacobian;
    colvec lyapunov;
    int p;
    int maxNumberOfPoints;
    double timeStep;
    double epsilon;
    //@}
};

#endif // DYN_LYAPUNOVTIMESERIES_H
